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Abstract 



A novel parallel algorithm for matrix multiplication is presented. The hyper-systolic 
algorithm makes use of a one-dimensional processor abstraction. The procedure 
can be implemented on all types of parallel systems. It can handle matrix-vector 
multiplications as well as transposed matrix products. 
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Matrix multiplication is a fundamental operation in most numerical linear 
algebra applications. Its efficient implementation on parallel high performance 
computers, together with the implementation of other basic linear algebra 
operations, is therefore an issue of prime importance when providing these 
systems with scientific software libraries [1]. Consequently, considerable effort 
has been devoted in the past to the development of efficient parallel matrix 
multiplication algorithms, and this will remain a task in the future as well. 

A general rule, which applies not only to matrix multiplication is, that the 
choice of a proper parallel algorithm strongly depends on the architecture 
of the parallel computer on which the algorithm is to run. System aspects, 
such as SIMD or MIMD mode of operation, distributed or shared memory 
organization, cache or memory bank structure, construction and latency of 
the communication network, processor performance and size of local memory, 
etc., may render an algorithm which is highly efficient for one system rather 
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impractical for another system. Even on a given system it may be necessary 
to switch algorithms in different problem size domains. 

As a consequence, one needs a diversity of algorithms for one and the same 
operation, as well as systematic design approaches which allow to construct 
new algorithms or to modify existing ones in such a way that they suit both 
a given implementation system and problem size domain. In the following, 
one such design approach is presented and applied to matrix multiplication. 
The procedure will lead us to a novel class of parallel matrix multiplication 
algorithms which are applicable to distributed memory computers whose inter- 
connection pattern includes a ring as a subset of the system connectivity. This 
novel scheme is called the hyper-systolic matrix multiplication, as it is based 
on the hyper-systolic parallel computing concept [2]. The latter is generaliz- 
able to any kind of commutative and associative operation on abstract data 
types [3]. The communication complexity of the hyper-systolic matrix multi- 
plication is 0(n 2 pz), with n being the matrix dimension and p the number of 
processors, and thus, it is comparable to best parallel standard methods. 

The work presented is part of a program to develop a novel practicable form 
of distributed BLAS-3 (PBLAS-3). In a forthcoming publication, we will give 
algorithms for transposed matrix products and level-2 linear algebra compu- 
tations. 

The paper is organized as follows: in Section 2, we shortly comment on systolic 
algorithms and the origin of hyper-systolic algorithms and review in brief 
the problem of distributed matrix multiplication. In Section 3, we develop 
the one-dimensional (ID) hyper-systolic matrix multiplication algorithm in a 
systematic way, starting from Fox' algorithm on a two-dimensional (2D) mesh 
[4]. In Section 4, the concept of the hyper-systolic algorithm involving two 
data arrays is introduced. The algorithmic presentation of the hyper-systolic 
matrix multiplication is given in Section 5, Section 6 deals with the mapping of 
the problem onto parallel systems and finally, Section 7 presents some results 
from implementations on a SIMD computer. 



2 Background 

2.1 Systolic arrays 

Systolic arrays are cellular automata models of parallel computing structures 
in which data processing and transfer are pipelined and the cells carry out 
functions of equal load between consecutive communication events. Systolic 
algorithms are parallel algorithms which, as far as abstract automata models 
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are concerned, make efficient use of systolic arrays. For more precise definitions 
of systolic algorithms and arrays and for many examples, the reader is referred 
to the monographs under references [5] and [6] (for a number of systolic ma- 
trix multiplication algorithms see Chapter 3 of [6]). The original motivation 
behind the systolic array concept was its suitability for VLSI implementation 
[7,8]. Only a few systolic algorithms, however, have been implemented in VLSI 
chips or hardware devices. With the advent of commercially available paral- 
lel computers, systolic algorithms have found an attractive implementation 
medium for they match the local regular interconnection structure as present 
in or easily implementable on many parallel architectures. 

Systolic algorithms can efficiently be implemented in the SIMD model, where — 
being synchronous and regular — they avoid time consuming synchronization 
operations. Apart from these implementation issues, a very attractive feature 
is the availability of methodologies for the systematic design of systolic al- 
gorithms. Projection of regular dependence graphs has evolved as one such 
technique [5,6,9-16]. 

As shown elsewhere [17-19], such structures can easily be transformed into 
data-parallel programs. The pattern of systolic arrays induces characteris- 
tic features into the respective data-parallel programs. In particular, a data- 
parallel program realizing a systolic algorithm consists of a sequence of iden- 
tical steps organized in a loop whose counter corresponds to the clock of the 
underlying systolic array automaton model. Further, the local regular intercon- 
nection pattern of a systolic array results in the use of only local synchronized 
communication in the respective data-parallel program as exemplified by the 
shift-type operations (e.g. cshift and eoshift). 



2.2 Hyper- systolic algorithms 

Hyper-systolic algorithms have been introduced for various numerical appli- 
cations in order to further reduce the communication expense of systolic al- 
gorithms. The advantage of the hyper-systolic over the systolic data flow first 
has been demonstrated for the case of so-called n 2 -problems, i.e., numerical 
problems that involve 0(n 2 ) computation events on pairs of elements in a 
system of n elements [2]. 

The systolic computation of n 2 -problems on a parallel computer equipped 
with p processors involves 0{np) communication events. The hyper-systolic 
algorithm can reduce the communication complexity to O^p^), as has been 
shown for a prototype n 2 -problem, the computation of all n 2 two-body forces 
for a system of n gravitatively interacting bodies [20]. This success makes us 
confident that hyper-systolic processing can be applied to a variety of numer- 
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ical problems which lead to n 2 computation events. An important application 
is found in astrophysics where the investigation of the dynamics and evolution 
of globular clusters is of prime importance [21]. Further examples of applica- 
tions are protein folding, polymer dynamics, polyelectrolytes, global and local 
all- nearest-neighbours problems, genome analysis, signal processing etc. [22]. 

Hyper-systolic algorithms are an extension of the systolic concept [2]. Similar 
to systolic algorithms, data processing and transfer are pipelined and the cells 
carry out functions of equal load. The three main differences are: (i) use of a 
changing interconnection pattern throughout the execution of the algorithm, 
(ii) use of multiple auxiliary data arrays for storage of intermediate results, 
and (Hi) the possible separation of communication and computation. The 
combination of these three features leads to reduction of the communication 
overhead. 

The use of a changing communication pattern is due to the communication of 
data by different strides along a ID ring in different stages of a hyper-systolic 
algorithm. The regularity of the communication pattern is retained, however. 
As an example of a regular but changing communication pattern one can think 
of an algorithm in which each processor of a ring communicates data to its 
first, second, fourth, etc., neighbours in the first, second, etc., steps of the 
algorithm, respectively. 

Auxiliary data arrays are needed for temporary storage of intermediate results. 
In conventional systolic algorithms such results are accumulated by shifting 
them from cell to cell. In the hyper-systolic algorithm they are generated and 
accumulated in place for many cycles, using multiple auxiliary data arrays, 
and are subsequently used to compute the final results. 

Use of a regular but changing communication pattern can be found in some 
(conventional) systolic algorithms, such as the systolic implementation [6] of 
Eklund's matrix transposition algorithm on a hyper-cube [23]. Use of multiple 
data arrays for solving specific tasks, e.g. problem partitioning, can also be 
found in systolic algorithm literature (see Chapter 12 in [23].), however, for 
purposes different of those aimed at in hyper-systolic algorithms. The unique 
combination of the two features mentioned above in the hyper-systolic concept 
aims at achieving new quality: substantial reduction of the communication 
overhead. 

2.3 Matrix multiplication on parallel systems 

In order to carry out matrix multiplications on distributed systems one has to 
take care for communication efficiency, parallelism and scalability of the im- 
plementation. For scalability of the implementation it appears very important 
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that the data layout chosen is preserved in the course of the computation with- 
out need for reordering of initial or result matrices. Furthermore, the data flow 
organization should not hinder the efficient usage of register, vector or cache 
facilities that largely dominates the overall performance. A third requirement 
is the general applicability in linear algebra tasks such as block factorization 
algorithms, where, additionally, matrix-vector multiplication must be carried 
out effectively^. It is a general observation that many parallel matrix mul- 
tiplication algorithms [24,25] meet the criteria of efficiency, parallelism and 
scalability only partially. For illustration, let us for instance consider one of 
the favorite matrix multiplication strategies, Cannon's algorithm [26], applied 
to n x n-matrix multiplication, A x B. 

Elements and of the matrices A and B are assigned to cells on a 2D 
grid, indexed by In a first step, Cannon's method involves a pre-skewing 
of A and B along the rows and columns, respectively. Data movement with 
different strides for each row/column is required, cf. Fig. 1. We note that this 
pre-skewing cannot be organized in form of a regular communication pattern. 
The second systolic step of Cannon's method consists of circular shifts of A 
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Fig. 1. Preskewing of A and B. 

along the rows and of B along the columns, followed by the computation event, 
cf. Fig. 2. We note that one can use blocks of elements as matrix entries. 

The merit of Cannon's algorithm lies in its memory efficiency, as it is possible 
to arrange the computation in such a way that no cell holds more than one 
element (block) of each matrix. Disadvantages of Cannon's algorithm are pre- 
skewing and the fact that the optimal data layout for matrix multiplication 
requires one-to-all broadcast operations applied to matrix-vector multiplica- 

1 Unfortunately it often happens that the matrix-to-processor mapping, as chosen 
for optimal BLAS-3 performance, fails in the case of BLAS-2 applications. 
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Fig. 2. Systolic phase of Cannon's algorithm. 

tion. Furthermore the layout of the result vector differs from that of the input 
vector [27]. 

3 Design of a ID Matrix Multiplication 



In this section, we develop a ID hyper-systolic matrix multiplication starting 
from a 2D algorithm that is related to the well known matrix product scheme of 
Fox [4]. In our systematic design approach the 2D scheme will be transformed 
to a ID systolic array representation. 

We have chosen a skew ordering as the fundamental representation of the 
matrices A, B and C. As a consequence, we will be able to carry out the com- 
putation fully in parallel without even the requirement for indexed addressing 
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functionality of the target parallel computer. Furthermore, no reordering is 
required in the course of the computation because the skew representation is 
preserved. 



3.1 Parallel Computational Problem 

Given the nxm matrix A and the mxn matrix B, the matrix-matrix product 
reads 



On a distributed memory machine, the elements of A and B have to be spread 
appropriately across the different processor nodes. 

In our design approach, we will first show how to multiply 4x4 matrices A 
and B on a 2D grid of nodes of size 4x4. Then, the algorithm is mapped onto 
a ID processing array with 4 nodes. 

In Section 5, we present the hyper-systolic algorithm for matrix multiplication 
of p x p matrices on an array of p processors. 

In Section 6, we consider the general case Eq. 1 of multiplying a nxm matrix 
A and am x n matrix B on a ID system of p processors. In order to map the 
systolic model onto a parallel implementation machine we choose hierarchy 
mapping with two different strategies, block and cyclic assignment. 



3.2 2D Matrix Multiplication 
3.2.1 Data alignment 

In the following we make use of the concept of abstract processor arrays (APA) 
as defined in HPF [28]. 

The grid of boxes shown in Fig. 3 is such a 2D APA on which the matrices 
A, B and C involved in the operation C = AB are aligned in a column-skewed 
fashion. 
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Fig. 3. Column-skewed distribution of the matrices A, B and C on a 2D APA. 
3.2.2 Semi- systolic algorithm 

Let the p x p matrices A, B and C be distributed on a p x p processor array 
and let for simplicity p — 4. 



The algorithm consists of p (4) steps as follows: in the first step, the matrix 
elements &i,i, &2,2j ■ ■ ■ ; &p,p of the matrix B in the first row of the APA are 
broadcast along the corresponding columns and subsequently are multiplied 
with the elements of the matrix A in the respective columns of the APA. The 
products (see Fig. 4a) are accumulated in the corresponding elements of the 
matrix C, according to the arrangement shown in Fig. 3. 

At the end of the step, the elements of the matrix A are circularly shifted by 
one position to the left along the rows and by one position downwards along 
the columns (compare the positions of the elements of matrix A in Fig. 4a and 
Fig. 4b). 

In the second step, the processors of the second row of the APA broadcast 
their corresponding elements of the matrix B in the respective columns of 
the APA, Fig. 4b. This operation is followed by the same multiplication and 
accumulation operations and circular shifting of the elements of the matrix A 
as in the first step. 



The algorithm proceeds with similar steps in which the processors of the third, 
through p-th row of the array broadcast in turn their elements of the matrix B, 
and the elements of matrix A are circularly shifted to the left and downwards 
by one position, cf. Fig.4c-d. After a total number of p such steps all partial 
products which belong to the elements of the matrix C are accumulated in the 
corresponding processors. 

The algorithm is called semi-systolic as it involves a broadcast operation of 
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Fig. 4. Computation of partial products on a 2D APA. 
the row elements of B in each step. 



3.2.3 Semi-hyper- systolic algorithm 

We next consider the semi-hyper-systolic variant which corresponds to the 
semi-systolic algorithm just described. The initial distribution of data is shown 
in Fig. 5. The distribution of the matrices A and C is the same as the one shown 
in Fig. 3. The distribution of the matrix B is obtained from the distribution 
shown in Fig. 3 by circular shifting of the elements of B in the i-th. row of 
the processor array by a stride of (i — 1) modx, with K = 2 for p = 4. In the 
particular example with p = 4, the second and the fourth row of B are shifted 
by one position. 

The algorithm consists of p (4) steps as shown in Fig. 6. In every step, each 



9 






Fig. 5. Initial distribution of data for the semi-hyper-systolic matrix product on a 
2D APA. 

processor multiplies the elements of the matrix A with the elements of the ma- 
trix B which it receives via the associated broadcast line. The partial product 
thus computed is accumulated in one of two local variables. The local vari- 
ables represent elements of two auxiliary arrays and which distributed 
across the processor array. They are used to accumulate partial results for the 
computations of the elements of the matrix C. Similar to the original systolic 
algorithm the processors in the first, second, . . . , p-th row broadcast the ele- 
ments of the matrix B they contain to all the processors of the corresponding 
columns in the first, second, . . . , p-th step of the algorithm, respectively, see 
Fig. 6. 

At the end of each step the matrix A is cyclically shifted downwards by one 
position. Unlike the original algorithm the horizontal shift of the matrix A is 
performed only every second step by a stride of two elements. 

The algorithm is completed by elemental addition of the auxiliary arrays 
and which is preceded by circular shifting of by one position to the 
left. 



3.2.4 Comparison to Cannon's algorithm 

The semi-hyper-systolic algorithm defined above can be compared with Can- 
non's algorithm. Each of the two algorithms has merits and shortcomings: 
Cannon's algorithm needs pre-skewing of the elements of the matrix A along 
the rows and the elements of matrix B along the columns of the 2D processor 
array which can be considered as a disadvantage since it gives rise to additional 
interprocessor communication that is not evenly distributed across the proces- 
sors. On the other hand, the operations to be executed by each processor in a 
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Fig. 6. The products shown in the upper and the lower halves of the processors are 
accumulated in the auxiliary arrays and C*- 2 -*. 

given step of the algorithm are rather simple, comprising one multiplication, 
one addition and two shift operations. Also, Cannon's algorithm needs only 
local processor interconnections. 

In contrast, the semi- hyper-systolic algorithm as illustrated in Fig. 6 needs 
broadcasting lines and more complex control of the operations which have 
to be carried out by the processors: in a given step, each processor has to 
broadcast an element to all other processors of the same column, and the 
matrix A is circularly shifted in two directions, where in horizontal direction a 
stride of 2 is required for every second step while in vertical direction a stride 
of 1 is carried out in each step. On the other hand, this algorithm avoids 
pre-skewing of the matrices. 



11 



3. 3 Matrix Multiplication on a ID Processor Array 

Next we transform the semi-systolic algorithm discussed above into a fully 
systolic one for a ID processor array, in which a given processor executes in 
sequence the tasks assigned so far to the processors working on a given column 
of the 2D processor array. This procedure will eliminate the disadvantages of 
the semi-hyper-systolic 2D array algorithm since no broadcasting lines are 
needed, and moreover, control structures become simpler. Downward shifts 
of A merely amount to a re-assignment within the systolic cell and do not 
involve actual communication between cells. Note that the application of such 
a transformation on the algorithm of Cannon will not have the same effect as 
it would still be required to pre-skew one of the matrices. 

3. 3. 1 Data alignment on a ID processor array 

In the case of a ID processor array, we assign each column of the 2D array, 
Fig. 3, to one processor. The resulting layout of the matrices A, B and C is 
shown in Fig. 7. The dotted lines indicate the location of the elements in the 
cells of the previous 2D array. 

a, 
c b " 

a 2 

r b 31 
°31 

b 41 

Wi 

Fig. 7. APA with ID systolic array alignment. 

3.3.2 ID systolic algorithm 

Again we start with the systolic computation. The algorithm needs p (4) steps. 
It does not require a broadcast of the matrix elements of the matrix B. In the 
first step, the matrix elements bi t i, 62,2, • • • , b PtP °f the m atrix B which reside 
in the first, second, . . . , p-th processor, respectively, are multiplied with the 
elements of the respective columns of the matrix A. The products (see Fig. 8a) 
are accumulated in the corresponding elements of the matrix C, according to 
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the arrangement shown in Fig. 7. At the end of the first step, the elements of 
the matrix A are circularly shifted by one position to the left along the rows. 
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Fig. 8. APA for the systolic computation of partial products on a !D array. 

In the second step, the second row of the matrix B is involved, see Fig. 8b. 
Its elements are multiplied by the elements of A and the partial products are 
accumulated in their proper locations, i.e., in the second step the products 
are copied to the elements of C which are circularly assigned one position 
downwards. Subsequently, A is circularly shifted to the left by one position. 

The algorithm proceeds with similar steps in which the elements of the third, 
through p-th row the matrix B are multiplied with the elements of A, the 
products being assigned to a row of C at locations i — 1 positions downwards 
the row in the i-th step, with the elements of matrix A being circularly shifted 
to the left, cf. Fig. 8c-d. After a total number of p (4) steps all partial prod- 
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ucts which belong to the elements of the matrix C are accumulated in the 
corresponding processors. 



3.3.3 ID hyper- systolic algorithm 

Let us turn now to the ID realization of the hyper-systolic algorithm. The 
initial distribution — as in the 2D case — requires a partial skewing of B along 
the rows, see Fig. 9. The distribution of the matrix B is obtained from the 
distribution shown in Fig. 7 by circular shifting of the elements of B in the 
i-th row of the processor array by a stride of (i — 1) mod^, with K = 2 for 
p — 4. 
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Fig. 9. Initial distribution of data for the hyper-systolic matrix product on a ID 
APA. 

Step by step each processor multiplies the element of the matrix A it contains 
with the corresponding element of the matrix B. The partial product thus 
computed is accumulated in one of two local variables [K = 2). The products 
computed in alternate steps are again stored alternately in one of the two local 
variables, in an analogous fashion as for the systolic algorithm, i. e., in the i- 
th step, the product is assigned to a row located i — 1 elements downwards, 
cf. Fig. 10. We emphasize again that A is shifted only every second step in 
horizontal direction to the left by two elements. For the general case with p 
processors, we have to carry a shift by K elements in every i^-th step. 

The algorithm is completed by elemental addition of the auxiliary arrays 
and C*- 2 -* which is preceded by circular shifting of by one position to the 
left. 
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Fig. 10. Hyper-systolic algorithm on a ID array. 
3.3.4 Summary 

In a systematic design approach, we developed a matrix multiplication algo- 
rithm on a ID abstract processor array, starting from an algorithm on a 2D 
APA. The algorithmic transformations led to a ID hyper-systolic scheme 

• that avoids broadcast lines required in the 2D case, 

• that, given p processors, shows a complexity of interprocessor communica- 
tion on the ID APA which is equal to that of Cannon's 2D algorithm, 

• that avoids skewing operations and reordering. 

So far, we illustrated the scheme using 4x4 matrices distributed column-wise 
over 4 systolic cells. In section 5, we generalize the 4x4 problem to a p x p 
system, distributed on a p processor array. In the case of the 4x4 matrices, the 
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shift constant was 2. For the general case, we shall introduce a stride K which 
is carried out K times, where the integer constants K and K fulfil KK = p. 



4 Hyper-Systolic Bases 

Before we present the general hyper-systolic matrix product, we define the 
hyper-systolic algorithm for two data streams. Such a situation is given in the 
computation of convolution and correlation problems. The scheme will later 
be adapted to the matrix multiplication [22]. 

4-1 Hyper- Systolic Recipe 

Let x and z be two ID arrays both of length n. Assume that functions 

F i = Qf(x i ,z j ) (2) 
j'=i 

for each i have to be computed, with © being an associative and commutative 
operator. The computation can readily be carried out by usage of systolic 
algorithms on a ring of systolic cells [20]. However, one can observe redundant 
interprocessor communication in this process that can be removed in hyper- 
systolic processing. The general recipe to find optimal hyper-systolic bases 
reads: 

Let the numerical problem be computable on two ID systolic processor arrays 
by use of a ID systolic algorithm. Let the two ID data streams x and z both 
of length n be mapped onto themselves by two sequences of (circular) shifts 
on the systolic array. The hyper-systolic scheme to compute the problem is 
constructed as follows: 

(1) For the systolic array x of length n, k replicas are generated by shifting 
the array x k times by strides a t , 1 < t < k, where all shifted arrays are 
stored as intermediate elements on the cells as arrays x t , 1 < t < k. 

For the systolic array z of length n, k' replicas are generated by shifting 
the array x k! times by strides b t >, 1 < t' < k', where all shifted arrays 
are stored as intermediate elements on the cells as arrays z t >, 1 < t' < k'. 

(2) The sequences of strides {a t }, 1 < t < k and {b t >}, 1 < t f < k' are 
determined such that 

(a) all pairings of data elements are present at least once, 

(b) the total communication cost is minimized. 
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(3) After each communication event the computations can be carried out and 
the results are assigned to the corresponding intermediate result arrays 
y t , 1 < i < & + 1. If elements occur more than once they are accounted 
for by a multiplicity table in order to avoid multiple counting. 

(4) One collector array y is moved by strides that follow the inverse of the 
sequence A k of strides of the initial phase. In each step of the back-shift 
phase the required intermediate result arrays y t are added to y. 

4-2 The Hyper- Systolic Optimization Problem 

Parallel machines support logical ID chains of processors in form of linear 
arrays or rings. However, circular shifts along the ID ring in general lead to 
different hardware communication expenses for different strides. 

The optimal sequence of strides for minimal interprocessor communication will 
depend on the interprocessor communication cost for a given stride. In order 
to minimize the communication cost effect on a given machine, we introduce 
a cost function C(oj), as a function of the stride a;. 

For the sake of argumentation, let us first assume the costs of communication 
for each array x and z on the systolic ring to be constant for any stride a iy hi. 
C(di) = C(bj) = const. 

Definition 1 - Optimization Problem for C(cii) = C(6j) = const. 
Let I be the set of integers m — {0, 1, 2 . . . , n — 1} G Nq, n G N. Find the two 
ordered multi-sets A k = (a® = 0,a±, 02, as, ... , a k ) G Nq +1 ofk + 1 integers and 
Bk' = (bo = 0, 61, 62, 63, . . . , bw) G N^ +1 of k' + 1 integers, with k + k' being 
a minimum, where each m G / , (0 < m < n — 1), can be represented at least 
once as the sum of two ordered partial sums 

m=(a i + Oi+i + . . . + a i+j ) + (b- + b l+1 + ... + b i+ .), (3) 

with 

0<i + j<k, i,je No, 0<i+j<k', i,j g N . (4) 



4.2.0.1 Lower bound on k + k' A lower bound for the minimal number 
of non-zero elements of A k can be derived that will deliver optimal complexity. 

Theorem 1 Let A k and be two bases solving the optimization problem for 
the hyper- systolic algorithm with 2 arrays. Then the minimal length k + k' is 
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given by 



k — k' — \pa — 1 . 



(5) 



PROOF. The total number of combinations required is n 2 as each element of 
the first array must come into contact with the n elements of the second array. 
Let the matrices 7i\ and 7^2 be realized by k — 1 and k' — 1 shifts respectively. 
In that case each element of the first matrix can be combined with k' elements 
of the second matrix, therefore the possible number of combinations will be 
nkk'. Given n = kk', the minimum number of circular shifts k + k' is attained 
for k = k! = y/n - 1. □ 



Therefore, the complexity for the interprocessor communication of a hyper- 
systolic algorithm for C(dj) = C(b,j) = const, is bounded from below by 3(y / n— 
1) shifts, where we have already included the costs for the back-shifts. 

4-3 C(oj) and C(6j) ^ const 

We now assume that the cost for a circular shift is a function of the strides a^. 
The optimization problem of definition 1 is modified only slightly, however, 
the construction of an optimal base can be quite complicated. 

Definition 2 - Optimization Problem for C(a,i) = C(6j) ^ const. 
Let I be the set of integers m = {0, 1, 2 . . . , n — 1} G Wq , n G N. Find the two 
ordered multi-sets = (ao = 0, oi, 02, 03, ... , a^) G Nq +1 of k + 1 integers 
and Bf = (b = 0, 61, b 2 , 63, ... , G of k' + 1 integers, with the total 

cost 



being a minimum, where each m G /, (t) < to < n — lj ; can &e represented at 
least once as the sum of two ordered partial sums 



k k> 



Ctotal = E C(aO + £ C(b-) 
i=i j =1 



i=l 



(6) 



m = + ai+i + . . . + a i+i ) + (ft- + b- +1 + ... + b- + -), 



(7) 



< i + j < k, 



i + j<k', !,jeN . 



(8) 
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4-4 Regular Bases 



The 4x4 problem presented in Section 3 uses so-called regular bases. This 
prescription turns out to be optimal for equal cost of any stride executed in 
circular shift operations on the ring. Regular hyper-systolic bases are advan- 
tageous as they require only two distinct strides. 

Definition 3 - Regular Bases The regular bases are given by 
A k=K ^: = (0,l, B kl=R „ v = (0,K,K,...,K) 

- v - ; - : A X A = 71.(9) 

K-l K-l 



The completeness of a base pair is defined in terms of the /i-range of the base, 
a notion borrowed from additive number theory [29,30]: 

Theorem 2 The h-range of a regular base is n. 



PROOF. Let 

r := mmod^ — > r < K, (10) 



as KK = n. There are K — 1 elements Oj = 1 G A k . Thus any r with 
< r < K — 1 r G No can be represented as partial sum by the elements 
a,i = 1, a« G .Afc. The partial sums of 5^', 

^^ = K^l = (j- ? + l)^<n, (11) 

i=j l=i 



are integer multiples of K. Adding the partial sums to r we can therefore 
represent any element m G /. Thus, the /i-range of the base pair Ak=K~i and 
B k , = x_ x is n, i. e., the base pair is complete. □ 

Theorem 3 The lower bound to the minimal length of the regular bases for 
a given h-range n is K = K = y/n. 

PROOF. The regular base A k is complete. 

k = K + K-l^k = K + — -l. (12) 
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Differentiation gives K = y/n. □ 



Theorem 4 The communication gain factor R that compares the regular hyper- 
systolic to the systolic algorithm is: 

R= n ~} (13) 
2K + K-Z 3 

PROOF. Let K = K = y/n. One needs K - 1 shifts by 1 and K - 1 shifts 
by K in forward direction and again K — 1 shifts by 1 in backward direction, 
respectively; therefore, the total number of shifts required turns out to be 

T=(2K + K-3). (14) 
The standard systolic computation requires n — \ shifts altogether. □ 



5 Hyper-Systolic Matrix Product 

Next we present the general formulation of the systolic and hyper-systolic 
matrix product in terms of a pseudo-code formulation. The size of the matrices 
is p x p and the ID processor array consists of p nodes. 

5.1 Systolic Algorithm 

The systolic version of the matrix product of two matrices A and B is given 
in Algorithm I. The matrices are represented in skew order. 

Algorithm 1 Systolic matrix-matrix multiplication. 

foreach processor i = 1 : p € systolic array 
for j = 1 : p 
for I = 1 : p 

ci,i = ci ti + cshift-colp~ J (a M ) b jti 
end for 
for k = 1 : p 

a k:i = cshift-row^(a fc)i ) 
end for 
end for 
end foreach 
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We have simplified the representation of data movements and assignments by 
introduction of two functions: 



cshift-row: horizontal circular shift of data by a stride of A; on a ring of cells 
numbered from 1 to n, cshift-row involves interprocessor communication. 



cshift-col: vertical circular shift by a stride of k for the vector of n elements 
within the systolic cells, cshift-col only amounts to memory assignments. 



The algorithm is completely regular. Each cell executes one compute operation 
together with an assignment followed by a circular shift of the matrix A in 
each systolic cycle. The skew order is not destroyed during execution of the 
algorithm. Note that for each processor inner cell assignment operations (col- 
cshift) are executed using equal strides in a given step of the parallel algorithm. 
Hence, one global address suffices and further address computations are not 
required! 



5.2 Hyper- systolic Algorithm 

It has been already noticed in section 3.3, that the complexity of the sys- 
tolic computation is not competitive with Cannon's algorithm. However, the 
hyper-systolic computation belongs to the same complexity class as Cannon's 
algorithm. It is fully pipelined and parallel, and does not require any skewing 
steps to align or re-align matrix elements. 



5.2.1 Regular Bases 

We employ the regular bases constructed for the hyper-systolic system. We 
add a third base C to account for the back shifts: 




(15) 




(16) 



k=K-l — 



(0, K, K, ...,K) 




(17) 



Ck'=K-l — 



(o, 1, 1, 
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5.2.2 Hyper- systolic matrix multiplication 

The hyper-systolic matrix multiplication, as given in Algorithm 2, proceeds 
within three steps. In the first part of the algorithm, matrix B is shifted K — 1 
times by strides of 1 along the systolic ring and stored as B l , < i < K — 1. 
However, as motivated above, for the case of matrix products, we can spare 
communication: it suffices to shift B in K row blocks of K rows each, where 
within each block the first row is shifted by a stride of and the last by a 
stride of K — 1. 

Algorithm 2 Hyper-systolic matrix multiplication. 

foreach processor i = 1 : p £ systolic array 

for j = 1 : p ! pre-shift of matrix B 

b j , i = c&fo-i*4 1 - i)modK \b j , i ) 

end for 

for j = 1 : K — 1 ! multiplication and shift of matrix A 

for I = 1 : K 
for n = 1 : p 

c n,i = cl n,i + cshift-colj, 1_0_1)X_il (an,i) b [{j _ 1)K+l]>i 
end for 
end for 
for I = 1 : p 

aij = cshift-row^(a M ) 
end for 
end for 
for I = 1 : K 
for n = 1 : p 

c n,i — c n,i + CSniTt COIp \a n ^) 

end for 
end for 

for j = 1 : K — 1 ! back shift and accumulation 

for I = 1 : p 

c M = c M + cshift-row p (c M ) 

end for 
end for 

end foreach 

After the preparatory shifts of B, the computation starts. K times, the multi- 
plication of A with K rows of the pre-shifted matrix B is carried out. After each 
step, A is moved to the left by a shift of stride K. The result is accumulated 
within K matrices C 
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Finally, the K intermediate result matrices C are shifted back according to 
base Cy while summed up to the final matrix C. The algorithm is very regular. 
The skew order is not destroyed during execution, and in any stage, only global 
addresses are required. 



5.2.3 Complexity 

The gain factor for the matrix product reads (note that matrix B is only 
partially shifted): 

Theorem 5 The gain factor R that compares the regular hyper- systolic matrix 
multiplication to the systolic algorithm is 

R= V -} «vg. (18) 
K+K-l 2 v 1 



PROOF. One needs 1 shift of the full matrix B, K — 1 shifts by K of matrix 
A and again K — 1 shifts by 1 of matrix C. Therefore, the total number of 
shifts required is 

T= (K + K-l). (19) 



The standard systolic computation requires p — 1 shifts of the matrix A. For 
the K = K = y/p, R « ^. □ 



5.2.4 Comparison to Cannon's algorithm 

In order to compare Cannon's algorithm and the hyper-systolic matrix multi- 
plication, we consider a 2D i/px ^processor array, where Cannon's algorithm 
is carried out, and a ID ring array of p processors where the hyper-systolic 
matrix product algorithm is computed. 

The total number of shift operations of Cannon's algorithm is 2^/p — 2 = 
2K — 2, while the number of shift operations for the hyper-systolic algorithm 
was 2K — 1. Thus, the complexities in terms of circular shift operations of both 
algorithms are equal for large p. We will see below how this fact translates 
into execution times on mesh and grid based machines. 
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6 Mapping on Parallel Systems 



So far we discussed the generic situation of the matrix dimension p being 
equal to the number of processors p. We now turn to the general case of 
square n x m-matrices with n,m > p. 

In order to map the systolic system onto the parallel implementation machine 
we choose hierarchy mapping of the systolic array onto the processors with the 
option for two different strategies, block and cyclic assignment. 



6. 1 Block Mapping 



The block assignment is applied in all standard algorithms, as it allows to 
exploit local BLAS-3 routines, like dgemm, by which a very high efficiency 
of local computations can be achieved. While a small block of the matrix A 
is hold in the cache or in the registers (thus avoiding cache-to-memory data 
transfer), in turn, only the columns of B and C must be exchanged, and all 
computations in which the given part of A is involved can be carried out. 
In this way, the ratio between the number of computations and the cache-to- 
memory traffic is minimized to nearly 

2/ 3 / 

(20) 



(3nn + n 2 ) 2 



floating point operations per word for real data, with I being the dimension of 
the sub-block. Asymptotically, the full speed of the CPU should be exploitable. 

A n x m-matrix M is divided into p x p blocks of size (j x j^j or (j£ x ^j, 
M^Mij, i=l,...,p, 3 = 1,..., p. (21) 



The multiplication of A and B proceeds via sub-matrix multiplication denoted 

as (<g>): 



k=l 
C = AB. 



Bfcj, i=l,...,p, j = l,...,p, 



(22) 



Altogether a system of p x p of such blocks is assigned to the p processor 
array. Now we can use each sub-matrix in the same manner as the scalar 
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matrix elements before. Therefore, the px p system of sub-matrices has to be 
row-skewed for A and column-skewed for B. 



6.2 Cyclic Mapping 



Each block is distributed across the processors as described above for the 
generic case. The blocking of the n x m-matrix M into x ^ blocks, leads 
to blocks of size p x p, 

M — > Mij, i 1 -. j 1 -. (23) 

p p 



The multiplication of A and B proceeds via block- multiplication, ((g)): 



*Ti p p 

C = AB. (24) 

A skew representation is required for all blocks separately. Cyclic mapping 
leads to a system of ^ x ^ systolic processes that run in parallel. 

Cyclic assignment allows us to reduce the memory overhead of hyper-systolic 
computations. In general, K full intermediate matrices C are necessary. Using 
cyclic mapping, one can organize the computation in such a way that only 
one row of the blocks of the intermediate matrix C must be stored in a given 
phase of the algorithm. All the required shifts of the given part of A can be 
carried out while this part of A will not be involved in a further computation. 
Eventually, the corresponding row of C is shifted back and accumulated. 

A second interesting feature of cyclic mapping is the possibility to distribute 
the very blocks. This approach is interesting for machines with a hybrid ar- 
chitecture like the proposed Italian PQE2000 system [31]. In this approach, 
the rows of A and the columns of B are assigned to different processors each. 



6.3 Block- Cyclic Mapping 



One can combine block and cyclic mapping in a hybrid scheme that combines 
the advantages of both approaches. A good strategy is to choose the block 
size of the block mapping such that it is optimal for "local" BLAS-3. For the 
cyclic part one ends up with blocks of size p x p, with the entries being the 
BLAS-3 blocks. 
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7 Implementation Issues 



The class of algorithm presented can be useful on any type of massively parallel 
system with distributed memory. Mesh and grid-based connectivities might 
benefit as well as large work-station clusters. 

In the previous section, we have presented complexities of interprocessor com- 
munication in terms of circular shifts, irrespective of the actual communica- 
tion time of a shift. This time is a constant on workstation clusters usually 
connected via Ethernet, therefore complexities in terms of circular shifts will 
translate into real time in a straightforward way. On mesh and grid based 
machines ID rings in general can be realized as a subset of their system con- 
nectivity. 

As an example taken from real life, we have implemented a level-3 PBLAS 
code on the APE100 parallel computer. So far, on APE100 lack of indexed 
addressing has hindered an effective scalable implementation of PBLAS [32]. 
In our approach, we made use of a combination of block and cyclic mapping. 
Thus, we are able to use local BLAS, exploiting the CPU with high efficiency, 
and to reduce the memory overhead of the hyper-systolic algorithm. 

On APE100, for real data, the optimal elementary blocks are of size 6x6. For 
complex data, the size is 4 x 4. The full matrix is blocked to p x p matrices 
which are distributed on the ring and the elements of which are the elementary 
blocks. Only the p x p matrices are skew, the elementary matrices remain in 
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Fig. 11. Performance of systolic and hyper-systolic PBLAS-3 (block-cyclic mapping) 
on a 128-node APE100, for real and complex data. 



26 



normal order. 



We have benchmarked a 128-node APE/lOOQuadrics QH1 using real and com- 
plex data. Fig. 11 shows the performance results. 

The theoretical peak performances (single node!) for Quadrics are 63 % for real 
data and 88 % for complex data, as can be inferred from the maximal ratio of 
computation vs. memory-to-register data transfer times. Hyper-systolic matrix 
multiplication leads to a peak performance of 65 % of peak speed, which 
translates into 75 % of the theoretical performance. 



8 Summary 



The ID hyper-systolic matrix multiplication algorithm is a promising alter- 
native to 2D matrix product algorithms. With equal complexity as standard 
methods, the hyper-systolic algorithm avoids non-regular communication and 
indexed local addressing. Hence, the hyper-systolic matrix product scheme is 
universal: it is applicable on any type of parallel system, even on machines 
that cannot compute local addresses. The method preserves the alignment of 
the matrices in the course of the computation. Besides the fact that trans- 
posed matrix products can be carried out on the same footing, alignments for 
the optimal hyper-systolic algorithm are efficient for matrix-vector operations 
as well. 

So far, in a feasibility study, the hyper-systolic matrix product has been imple- 
mented on APEIOO/Quadrics systems [33,34]. We will present details of the 
implementation on this system for matrix and transposed matrix products in 
Ref. [35]. 
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